TRANSPORT IN TRANSITORY DYNAMICAL SYSTEMS * 

B. A. MOSOVSKY + AND J. D. MEISS + 

Abstract. We introduce the concept of a "transitory" dynamical system — one whose time- 
dependence is confined to a compact interval — and show how to quantify transport between two- 
dimensional Lagrangian coherent structures for the Hamiltonian case. This requires knowing only 
the "action" of relevant heteroclinic orbits at the intersection of invariant manifolds of "forward" and 
"backward" hyperbolic orbits. These manifolds can be easily computed by leveraging the autonomous 
nature of the vector fields on either side of the time-dependent transition. As illustrative examples 
we consider a two-dimensional fluid flow in a rotating double-gyre configuration and a simple one- 
i—t and-a-half degree of freedom model of a resonant particle accelerator. We compare our results to 

f**") those obtained using finite-time Lyapunov exponents and to adiabatic theory, discussing the benefits 

^vj and limitations of each method. 
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1. Transitory Systems. Invariant manifolds have long been recognized as im- 
portant structures that govern global behavior in dynamical systems. Hyperbolic 
manifolds in particular, by their very definition, relate information about the ex- 
ponential contraction and expansion of nearby trajectories within the flow and so 
play crucial roles in the dynamics of such systems, lending insight into the mecha- 
nisms by which chaos, mixing, transport, and other complex global phenomena occur. 
Transverse intersections of stable and unstable manifolds give rise to lobes defining 
of packets of trajectories that exit or enter coherent structures or resonance zones 

^i bounded by pieces of the manifolds. Thus, tracking these lobes provides a means for 

h> quantifying flux between coherent structures in the flow. 

\& The treatment of mixing and transport for aperiodically time-dependent flows, 

however, requires the development of new methods because the concept of invariance 
may be too strong and may not even lead to physically relevant structures. One 
popular method, the finite-time Lyapunov exponent (FTLE) has been used extensively 
in recent years to compute local approximations of invariant manifolds and identify 
structures that remain coherent in the Lagrangian sense on some finite time interval 
pTl 152] [30 ] . Another idea uses a nonautonomous analog of hyperbolic orbits, called 
a distinguished hyperbolic trajectory [HI HZl 1231 HI] j and a third technique identifies 
approximately invariant regions as eigenfunctions of the Perron-Frobenius operator 



While these Lagrangian frameworks for identifying coherent structures in nonsta- 
tionary systems have seen broad application, using them to accurately quantify trans- 
port and mixing over finite timescales remains a difficult task. A few recent studies 
have managed to give numerical estimates for finite-time transport and mixing within 
aperiodic time-dependent flows [IT] QTJ [51 03]. In addition, the instantaneous flux 
across a "gate surface" can also be estimated using only local Eulerian information 
pm Hj , and the results have been applied to geophysical flows defined both analyti- 
cally and by discrete data sets. One such application is the investigation of eddy-jet 
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interactions, in which the strength of the interaction is quantified by the amount of 
fluid entrained by or ejected from the eddy [3H1 US]- However, the accuracy of these 
methods is tied to the rate of change of the Eulerian velocity field, and, moreover, the 
instantaneous flux through a surface does not provide an estimate for the finite-time 
transport between two disjoint coherent structures. 

Since the mere identification of Lagrangian coherent structures (LCS) in fully 
aperiodically time-dependent systems is a challenging problem itself, we study here a 
simpler problem — the quantification of finite-time transport between coherent struc- 
tures in two-dimensional systems undergoing a transition between two steady states. 
Our methods are Lagrangian, in the sense that they rely on knowing certain key tra- 
jectories, and as such they allow us to compute the transport between two bounded 
coherent structures within the flow, as opposed to knowing only the flux across a sin- 
gle gate surface. Even though we restrict the time-dependence to a compact interval, 
this problem provides, we believe, some insight into the more general aperiodic case. 

The systems of interest to us here are stationary in the limits t — > ±00; these were 
called asymptotically autonomous systems by Markus [41 . However, we will assume 
that the system is nonstationary only on a compact interval, and will refer to these 
systems as transitory: 

Definition 1.1 (Transitory Dynamical System). A transitory dynamical system 
of transition time r is one that is autonomous except on some compact interval of 
length t, say [T p ,Tf] with Tf — T p = r. Thus on a phase space M, a transitory ODE 
has the form 

x = v ix ,t), v {x ,t) = [ p F \i] <<;*, (i.i) 

where P : M — > TM is the past vector field, F : M — > TM is the future vector field, 
and V(x,t) is otherwise arbitrary on the transition interval [t p ,t/]. 

Though P and F are assumed autonomous, they could have arbitrarily compli- 
cated dynamics. The transition time t = Tf — t p is, relative to the time scales of P 



and F, an especially important parameter for ( 1.1 ) and, without loss of generality, we 
may set t p = so that Tf =r. 

It is not hard to envision many physical situations to which transitory dynamics 
would pertain. For example, any dynamical system that depends upon some param- 
eters, x = V(x;p), can be made transitory if the parameters become time-dependent, 
p — > p(t), and are allowed to switch from one state to another over a time interval of 
length t. We will consider a simple model of a particle accelerator in §3.2| for which 
this is the case, but more realistic models could also be used [Hj. A similar phys- 
ical system — also Hamiltonian — corresponds to a point particle in a billiard whose 
boundary evolves over time; for example, an ellipse with eccentricity that changes 
over a compact interval from one value to another (periodically oscillating boundaries 
have been studied in a number of cases, e.g. [31]). Ecological models could also be 
transitory if the environment undergoes a shift (e.g., leading to a change in carrying 
capacity); one could study the "transport" between basins of different equilibria under 
such a shift. Another class of examples corresponds to the change in flow regimes for 
a fluid in which there is an instability that grows and saturates. This could be due to 
external forcing (e.g., Rayleigh-Bernard convection or Taylor-Couette flow), to a flow 
that is driven by chemical reactions (e.g., Marangoni flow), or to an instability leading 
to eddy creation (e.g., for quasigeostrophic flows 47J). Finally, many models of fluid 
mixing could be studied in transitory regimes; for example, laminar flow through a 
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pipe with a finite number of bends between two straight sections could be considered 
transitory [27J 13] and a fluid-filled cavity [TU1 [32J [2] could be driven with a transitory 
mixing protocol. 

For more general aperiodically time-dependent systems, infinite-time information 
cannot be used to identify coherent structures: indeed, in many fluid mechanical and 
observational applications, the behavior of the system is known only on a finite in- 
terval. This requires alternate definitions of approximate invariant manifolds. These 
can be defined by extending the vector field to an infinite-time domain; however, 
since the extension is not unique, neither are the resulting manifolds [201 ED [HI \EM ■ 
Nevertheless, if the time of definition is long compared to the local expansion rates, 
this nonuniqueness results only in exponentially small corrections. By contrast, for 
transitory systems, stationarity outside the transition interval [0, r] implies the classi- 
cal notions of stable and unstable invariant manifolds can be used — see below. Thus, 
transitory systems provide an aperiodic, time-dependent setting for which unique in- 
variant manifolds exist. The advantages of using this asymptotic information will 



become clear when contrasting our methods with those employing FTLE in £3.1 and 
§3.2[ In any case, a more general vector field that is defined only on a finite interval, 
[0,t], say, could be extended by adding autonomous past and future vector fields to 
give a transitory vector field. Thus the relative unimportance of the form of the vector 
field's extension, given a long enough interval of definition, also applies to our case. 

While we are not aware of other studies of transport for transitory systems, some 
aspects of transport for asymptotically autonomous systems in which the forward and 
backward limits are identical, 

lim V(x,t) = G(x), (1.2) 

t— ^±oo 

have been considered by Wiggins and collaborators j40j [50] . A well-studied case 
corresponds to the adiabatic limit, when r is large compared to the dynamical time 
scales of P and F [25]. Indeed, when < |1.1| ) is Hamiltonian and adiabatic, then the 
actions of the frozen system P (or F) should be approximately preserved. However 
as is well-known, adiabatic invariance breaks down near separatrices of the frozen 
system [5J [45] and cases in which these separatrices sweep through a large portion of 
the phase space during the transition are of most interest to us. As was first noted 
by Elskens and Escande [15) . when the time-dependence is periodic the entire region 
swept by the separatrix becomes a "lobe" in the adiabatic limit. Transport properties 
have been studied in this limit 25, 5, 26 ; however, these researchers assume that the 
frozen system has a parameter dependent curve of hyperbolic equilibria and this is 



typically not true for ( 1.1 ). Moreover, we will not assume that r is large. 

In this paper we are interested in quantifying the transport between coherent 
structures of the past and future vector fields P and F for two-dimensional systems 



of the form (1.1 1. To define these, it is natural to consider hyperbolic orbits of P 
and F and their stable and unstable manifolds since initial segments of these often 
define invariant or nearly invariant structures such as resonance zones [391 113| 134] . 
However, determining which structures of P and F are relevant to the dynamics of 
the full nonautonomous vector field V requires special attention. 



It is natural to think of the dynamics of ( 1.1 ) as occurring on the extended phase 
space Mxl. We will assume that V has a complete flow, (fit 1 ,t '■ M — > M for any 
to,t\ € R, where <fti,t rnaps a point from its position at t = to to its position at 
t = t\. Then every point (x, t) £ M X K has an orbit 

7(z,to) = {(Vti,* (»)»*i) : ^i e K} C M x R, 
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and such sets are invariant in the sense that for any r£l, j{x, t) = 7((y9 r-t (a;), r). Of 
course, this notion of invariance is not restrictive: any subset S C M can be taken 
to be the time-r slice through an invariant set j(S,t) — {(x,t) : x £ if t . T (S),t £ R}. 
More generally the time-i slice of any set TV £ MxR can be defined using the standard 
projection 7r : M x R — >• M by 

N t = Tr(ND{{x,t):x£M})cM. (1.3) 

Thus 74 (S,t) = S. 

If A is any invariant set of the past vector field P, then the set {(A, t) : t < 
0} C M x R is a backward-invariant set in the extended phase space. For example, 
equilibria and periodic orbits of P are slices of backward-invariant sets of V. Similarly 
any invariant set of F is a time-i slice of a forward-invariant set of V for each t > r. 



Since the nonautonomous portion of the dynamics of (1.1 1 is assumed to occur on 
a compact interval, it can be effected by a map, the transition map T : M — » M, 
defined as 

T(a;) - ^ T ,o(x). (1.4) 

Consequently an invariant set A of P becomes T(A) at time r and thereafter evolves 
under F. If the dynamics of P and F are known, then the only nontrivial work we 
must do is to characterize the map T. 

In addition to equilibria and periodic orbits, the unstable manifolds of orbits of P 
and stable manifolds of orbits of F are also slices of invariant sets for V. For example, 
let W a (A, X) denote the unstable manifold of an invariant set A under a vector field 
X; it is the set of points that are asymptotic to A as t — >• — oo. Consequently, if A is 
an invariant set of P, its unstable manifold is a slice of the unstable manifold of the 
invariant set 7(A, 0) of V: 

W"(A,P) = IF t u ( 7 (A,0)) = T4?( 7 (A,0),n t < 0[] 

However, a stable manifold W S (A, P) is not a slice of a stable manifold for V since, 
due to the transition, it does not in general consist of points asymptotic to the orbit 
of A as t — > oo. Instead, since the forward dynamics for any t > r is determined 
by F(x), the stable manifolds of invariant sets of F are time-i slices of the stable 
manifolds for V. The unstable manifolds of invariant sets of F have little dynamical 
relevance to V. 

An example is sketched in Fig. |1.1| Here p is a hyperbolic saddle for P, but its 
orbit under V, j(p, 0), is not hyperbolic. Indeed, when the point T(p) does not lie 
on a hyperbolic orbit of F, -f(p, 0) has no stable manifold. Nevertheless, it does have 
an unstable manifold W u (7(^,0)) and the temporal slices of this manifold coincide 
with W u (p, P) when t < 0. Furthermore, this manifold is an invariant set of V 
with W^(^{p : 0)) = T(W n (p, P)) and its subsequent structure is obtained by simply 
evolving this set with F. We will say that such an orbit is backward hyperbolic. 
Similarly, if / is a hyperbolic fixed point for F it is a forward hyperbolic orbit of V, 
but not generally hyperbolic under V. Slices of the stable manifold W s (j(f, r)) agree 
with W s (f, F) for each t > t. We formalize these notions as follows: 

Definition 1.2. A time-to slice of an invariant set A of a transitory dynamical 
system is backward hyperbolic provided <Pt,t (At ) * s hyperbolic under P for all t < 0. 
It is forward hyperbolic provided yt,t (At ) is hyperbolic under F for all t > r. 



1 As indicated here, we often will omit the V from the notation as it represents the default 
evolution. 
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Fig. 1.1. A backward-hyperbolic orbit ~f(p, 0) with homoclinic loop T p C W u (p,P) and a 
forward-hyperbolic orbit "f{f, r) with homoclinic loop Tf C W s (f,F). Here the unstable manifold 
of the orbit of p intersects the stable manifold of the orbit of f at t = r at the heteroclinic points 
hi, h'2 <Z T(T p ) fl r^ under a transitory flow. 



In the simplest case, dim(M) = 2 and T p c W u {p,P) and T f C W s (f,F) are 
homoclinic loops for P and F, respectively, as sketched in Fig. |1.1| The regions 
bounded by T p and Tf can be thought of as Lagrangian coherent structures (LCS). If 
the former evolves under the map ( 1.4 ) to intersect the latter, then we can characterize 



the transport from one coherent structure to the other by these intersections. 

Lagrangian coherent structures for nonautonomous systems are usually defined 
in terms of finite time stability exponents. Haller and Yuan define them as regions 
bounded by "material lines with locally the longest or shortest stability or instability 
time" [21], while other researchers refer to LCS as curves or surfaces instead of regions, 
and identify these as "ridges" in the finite-time Lyapunov exponent field 52 , 42 , 6 . In 
either case, since these coherent structures arc defined only by finite time information, 
they may also persist only for some finite time. For (1.1 1, we will think of LCS as 
being objects bounded by separatrices of P for t < or of F for t > r. These 
structures may also be ephemeral under the vector field V; however, for a transitory 
system we think of only one event, encapsulated by the transition map T, as creating 
or destroying an LCS. 

In Fig. |1.1[ the coherent structure bounded by T p in the past vector field P is 
destroyed by the transition map T, giving rise to a new structure bounded by Tf in 
the future vector field F. There are two heteroclinic points {hi, /12} = T(T P ) n Tf 
in the time-r slice; they are backward asymptotic to p, forward asymptotic to /, and 
hence fully hyperbolic under V. Consequently, the set of orbits that begin inside T p 
but evolve to escape from Tf is defined by the lobe R bounded by the segments of Tf 
and T(T p ) between hi and hi- More generally, there may be more than one lobe and 
each lobe boundary may consist of more than two manifold segments, or equivalently, 
contain more than two heteroclinic points. We discuss the computation of lobe areas 
in this general case in §2.1| 

For higher dimensional flows, even though the invariant manifolds W s {p, P) and 
W u (f, F) are not dynamically relevant for the full vector field V(x, t), they may useful 
for defining resonance zones of P and F, where a resonance zone with a "small" 
escaping flux is a "nearly" invariant set [SH] H31 GH] ■ Particles initially in a resonance 
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zone of P when t < will be approximately trapped up to time 0, and those that find 
themselves in a resonance zone of F at t = r will be approximately trapped in the 
future. This makes low-flux resonance zones good candidates for coherent structures 
in the past or future vector fields. 

In the remainder of the paper, we consider several simple examples of transitory 
systems for which the full vector field is a convex combination of the past and future 
vector fields: 

V(x, t) = (1 - s(t))P(x) + s(t)F(x). (1.5) 

Here s : R — > [0, 1] is a transition function satisfying 



for transition time r. While ( 1.5 ) is transitory in the sense of Dcf. 1 1 . 1 1 for any function 



s satisfying (1.6), for simplicity we usually take s to be monotone nondecreasing. 
Examples of such functions of varying smoothness are given in App. [K\ 

As a first example, consider the traveling wave model studied by Knobloch and 
Weiss [28[ 153] . This model for an incompressible, two-dimensional fluid is given in 
terms of the stream function 

ip(%,y,t) = i>o(x,y) +eb(t)ipi(x,y), 

tpo(x,y) = -cy + Asm(kx)s'm(y), (1.7) 

i>i{x,y) = y, 

on the domain M = [0, 2ir] x [0, tt]. The vector field is given by V — z x VV>, where z 
is the unit normal to the ccy-plane, so the equations of motion are 

± = -J^ v = §-J, (i.8) 

This system is Hamiltonian with (x, y) representing the coordinate and momentum 
and H(x,y,t) = -ip(x,y,t). 

While Knobloch and Weiss studied the dynamics of -00 with time-periodic per- 



turbations, an asymptotically autonomous case of (1.7) was studied in [40j [50]. The 
latter assumed that lim t ^ ±00 b(t) — > b^ € E so that the past and future vector fields 
are equal. For this case b(t) is a "bump function" that plays the role of the transition 



function s in (1.5). If instead, we assume that b{t) has support only on the compact 
interval [0,r], then (1.8) is transitory in the sense of Def . [P] One such function is 



b(t) = s(t)(l-s(t)), (1.9) 



where s(t) is any transition function (1.6). 

For ( |1.7| ) with (1.9), P is identical to F, and when \c\ < \A\ there are two hyper- 
bolic equilibria on each of the lines y = and y — tt (see Fig. |1.2[ ) . The four saddles 
of the past and future vector fields are denoted p.; and fi, respectively. Note that 
as subsets of M, pi = fa however, as subsets of the extended phase space, they are 
points on different temporal slices and so their evolution under V is distinct. For the 
simple perturbation ifri , the orbits of these points under V remain on their respective 
horizontal lines. Moreover, if b(t) > 0, then <fit,o(pi,2) — > /i and ft,o(P3,i) — > /3 as 
t — > oo. 
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Fig. 1.2. Frozen time stream function at t = (left) and t = T (right), and coherent structures 
for | fj.7| ) with A = k = e = 1 and c = 0.5. The slices of the unstable manifolds in the right pane are 
obtained by numerically computing the transition map for r = 4 using the bump function (1.9ty with 
s(t) = s\(t/r), the cubic transition function in \A . w . 



The vector fields P and F each have two resonance zones in M and these form the 



coherent structures of interest for the system (1.7 1. For P, the two resonance zones 



should be thought of as being bounded by unstable manifolds since these will be slices 
of unstable manifolds of V; one is bounded by branches of W u (jpi, P)UW u (p2,P) and 
the other by W u {p%, P) U W u (p4, P), as shown in Fig. |1.2| For F the resonance zones 
are bounded by stable manifolds: W s (fi,F) U W s (f 2 , F) and W s (f 3 , F) U W s (fi, F). 
The unstable manifolds at t = evolve according to (1.8), and for the case shown in 



Fig. |1.2| numerical integration indicates that these intersect the stable manifolds at 
two heteroclinic orbits, 



fci(t) = W?(pi) n W/(/ 2 ), 

h 2 (t) = w t u ( P3 )nw t s (n). 



(1.10) 



The lobes formed by these intersections, labeled R\ and R 2 , correspond to the trajec- 
tories that begin inside the resonance zones of P and end outside the resonance zones 
of F. By calculating the areas of these lobes we can quantify transport into and out 
of the coherent structures in the phase space, giving insight into the global dynamics 
of the system. 

2. Transitory Flux: Hamiltonian Case. Coherent structures, by their very 
definition, denote regions of the phase space in which nearby trajectories behave sim- 
ilarly; examples include vortices or recirculation regions in fluid flows and resonance 
zones in the phase space of Hamiltonian systems. Since these coherent structures 
are typically composed of a large number of trajectories and since their boundaries 
separate dynamically distinct regions in the phase space, they can provide quite a 
bit of information about the global behavior of the system. In particular, knowing 
the incoming and exiting flux helps paint a global picture of the Lagrangian effects 
of time dependence within the vector field. For example, for a particle accelerator we 
may wish to quantify the phase space volume corresponding to stable acceleration of 



particles (see [3.2). In this case, computing the flux between coherent structures of 
the pre- and post-acceleration vector fields gives the desired quantity. 

We proceed in this section to derive formulas for the flux between regions of 
phase space bounded by invariant manifolds of nonautonomous, one degree-of-freedom 
Hamiltonian systems (i.e. a system with 1~ degrees of freedom). The formulas ob- 
tained are the generalizations to the nonautonomous case of the action-flux formulas 
of [35] • For autonomous flows, these formulas were first obtained in [3S1[371|3S] and 



8 B. A. MOSOVSKY AND J. D. MEISS 

the nonautonomous case was studied in [25] in the adiabatic limit. 

We note that [50, 40 do compute lobe areas for asymptotically autonomous vector 



fields, in the sense of ( 1.2 ). However their theory relies on P = F and that these vector 
fields have a saddle equilibrium with a homoclinic trajectory. In the theory we present 
here, the saddles of P and F need not be the same and there need not be a homoclinic 
orbit of either autonomous vector field. Instead we focus our attention on heteroclinic 
orbits of the time-dependent vector field V. 

2.1. Flux by the Lagrangian Action. Here we will compute the flux for a l| 
degree-of-freedom Hamiltonian vector field V — (d y H, —d x H) for a nonautonomous 
Hamiltonian H(x, y, t). More formally, if u is a symplectic formF] e.g. u = dx A dy, 
then the Hamiltonian vector field is determined by zyu = dH . In this section we do 
not need to assume that the system is transitory, but we do assume that the phase 
space M is two-dimensional and that the symplectic form is exact. 

As discussed in |T] computing the flux corresponds to finding the area of some 
closed and bounded region R C M. Using Stokes' theorem, the resulting two- 
dimensional integral can be immediately reduced to an integral over the boundary: 

Area(i?) = / u = - / u, (2.1) 

JR JdR 

where v is the Liouville one-form defined by u = — dv; for example, v = ydx. When 
R is bounded by segments of stable and unstable manifolds, the flux formulas of [35] 



reduce integrals of the form (2.1| to action differences between orbits lying at the 
endpoints of the manifold segments. The action is given by an integral of the phase 
space Lagrangian, L : M x R — > R 

L(x,y,t) = tvv-H(x,y t t), (2.2) 

or, L = yx — H with x(x, y, t) = d y H{x, y, t). 



The simplest case is sketched in the left pane of Fig. 2.1 Suppose that 7/ is a 
forward hyperbolic orbit, j p is a backward hyperbolic orbit and R is a region in the 
time-r slice that is bounded by a stable-unstable pair of segments of time-r slices, 
S C W®(7/) and U C W^(j p ), that intersect only on their boundaries, ho and hi. 
Choosing the orientation of U and S consistent with a counterclockwise traversal of 



the boundary of R gives dR =U + S and dS = —dU = ho — hi. Then (2.1 ) yields 



Area(i?) = - I / v + / v\. (2.3) 

Therefore, to find the area we may compute the integral of the Liouville form along 
segments of stable and unstable manifolds. The following two lemmas give formulas 



for computing these integrals in terms of the phase space Lagrangian (2.2 1 and the 
heteroclinic orbits hi(t) = tp tT (hi). 

Lemma 2.1. Suppose that the orbits of ho and h\ are backward hyperbolic and 
backward asymptotic, and that hi C W^(^f(hi,T)) is the time-r slice of the unstable 
manifold that connects these points, with dhi = hi — ho- Then 

[ p = AA-(ho,h x )= f [L(hi(s),s)-L(h (s),s)]ds. (2.4) 

JU J-00 



Our notation is given in App. B 
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Fig. 2.1. Sketches of time-r slices of lobes formed by stable and unstable manifold segments 
S and U. In the left pane the lobe is bounded by a pair of segments, and in the right pane by two 
pairs. Each segment is bounded by a pair of bi-asymptotic heteroclinic points hi and hi+i. 



Proof. If we differentiate v along V, Cartan's homotopy formula (B.7| and (2.2| 



gives 



dt 



V = CyV 



-lyW + d{lyV) = dL, 



where Cy is the Lie derivative (B.5|. Integrating this from t to r, using (B.8| gives, 
for any t, 



v-Vt,r v = 



t v = I foV**,r v d s = / ^Ir^ds. 



Consequently 



(J d(<p* SjT L)\ ds + J v * tT u 



dL ds ■ 
[L(hi(s),s)-L(h Q (s),s)]ds + 



(2.5) 



'Vt,r(W) 

Since ho(t) —t hi(t), the length \(p t _ T (U)\ — > as t — >• —00. Taking this limit yields 



the result (2.4 1. D 



We note that AA~(ho,hi) in (2.4| is the difference between the "past actions" 
of the orbits of ho and hi. A similar result, with an important sign change, holds for 
stable segments S C W*(~f(f, r)) connecting hi to h a . In this case we replace U by S 
and let t —¥ +00 in (2.5 1 to obtain 



Lemma 2.2. Suppose that the orbits of hg and hi are forward hyperbolic and for- 
ward asymptotic, and that S C W^( , y(hi,T)) is the time-r slice of the stable manifold 
that connects these points, with OS = ho — hi . Then 



v = -AA+(hi,ho) 



[L(h (s),s) -L(hi(s),s)]ds. 



(2.6) 
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Here AA^(hi,ho) is the difference between the "future actions" of the orbits of 
ho and hi . The algebraic area of the lobe R in the left pane Fig. |2.1| can be obtained 



by plugging the results (2.4) and (2.6) into (2.3), yielding the action difference, 



/oo 
[L(h (s),s) - L(hi(s),s)]ds. 
-OO 



(2.7) 



Note that r appears nowhere in this result; indeed, the lobe area is independent of 
the time at which it is measured since the Hamiltonian flow is area-preserving. 



It is important to note that equation ( 2.4 ) can be used to compute the area under 



any segment of an unstable manifold between two backward-asymptotic points ho and 



hi. That is, the points need not be bi-asymptotic as in Fig. 2.1 For example, if h is 



replaced with p(r), equation (2.4) can be directly used to integrate the Liouville form 
v along the initial segment of W r "(7 P ) from p(r) to hi. Similar generality holds for 
equation (2.6). Moreover, the derivations of (12. 41) and (2.6) do not require that the 



system be transitory; they are valid for any nonautonomous Hamiltonian vector held. 
Lobes can have a more complicated structure, even in the autonomous case [48] ; in 
general, a lobe at time r is a region R bounded by an alternating sequence of stable 
segments of 7/ and unstable segments of 7 P whose intersections are topologically 
transverser] An example with two pairs of segments is shown in the right pane of 
Fig. |2.1l The general formula for the area of such a lobe follows easily from ( 2.4 ) and 



(2.6). Suppose that R is bounded by 2N such segments, and label the heteroclinic 
points hi, i = 0, . . . , 2N — 1 in a counterclockwise ordering on dR setting h2N = ho- 
Without loss of generality, suppose the segment joining ho and hi is a portion of 
unstable manifold, call it U . Again using a counterclockwise ordering, label the 
next segment So, followed by Ui, etc., so that dR consists of the alternating sum 
of N unstable and stable segments Ui and Si, i = 0, . . . , N — 1. The orientation of 
these segments is chosen to be consistent with the counterclockwise boundary so that 
dUi — /i2i+i — h%i and dSi = h-n+2 — foi+i- Note that these orderings, as sketched 
in Fig. |2.1| are not necessarily the same as ordering along the manifolds W u or W s . 



Using (12. 41) and (2.6) for each of the integrals along dR, gives the lobe area 



N 



N-l 



Aiea(R) =^AA(^i,M = - ^ AA(h 2i ,h 2l+ i) 



(2.8) 



where AA is defined by (2.7). While the two sums in (2.8) are trivially identical, the 



first can be thought of as the sum of the action differences between orbits bounding 
segments of stable manifold along dR and the second as the negative sum of action 
differences between orbits bounding segments of unstable manifold. 

Lemmas |2.1| and |2.2| apply to arbitrary nonautonomous systems provided only that 
the points bounding the segments U and S are past or future asymptotic, respectively. 



For the special case of transitory systems, the integrals in (2.4) and <j2.6|) can be further 



simplified since the phase space Lagrangian is autonomous outside (0, r) : 



L{x,y,t) = 



L P (x,y) t<0 
L F (x,y) t>r 



3 If a heteroclinic point h arises at an intersection that is not topologically transverse, the two 
boundary segments adjacent to h can be combined (they necessarily have the same stability type) 
and h can thus be ignored when calculating the lobe area. 
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For example consider the lobe depicted in the left pane of Fig. 2.1 so that ho(0) and 
hi(0) both lie on W u (p, P), the unstable manifold of a hyperbolic fixed point p = p(Q) 
of P. Since this manifold is stationary for all t < 0, these points are merely time-shifts 
of one another under the flow of P. If we suppose that hi(Q) is further along W u (p, P) 
from p than ho(0) (as in Fig. 2.1 ), then there exists a tp < such that hi(tp) = ho(0) 
and thus 

hi(tp + a) = ho(a) Va e (— oo,0]. 

Consequently the integral in (|2.4|) reduces to two integrals over compact intervals: 



I v = AA-(h ,h l )= f L P (hi(s))ds+ f [L(hi(s),s)-L(h (s),s)]ds. (2.9) 

JU Jt P Jo 

Similarly, assuming that points ho and hi are oriented along the stable manifold as 
in Fig. 2.1 there is a tp > r such that h (tp) = /ii(t). Then (2.6| reduces to 

= -AA+(hi,ho) = - / r L F (h (s))ds. (2.10) 



Combining these yields the simplified form of the area (2.7) for the case of transitory 
systems: 

Area(i?) = / L(h (s) , s)ds - f Lp(h x {s))ds. (2.11) 



The formulas (2.9 1 and (2.101 prove useful in numerical computations, and they can 
be applied to each action difference of (2.8 1 in the case of a more complicated lobe 
bounded by 27V manifold segments. We will use them in the examples of S]3] 

2.2. Flux in the Adiabatic Limit. As in Def. |1.1| a Hamiltonian H(x,y,t) is 
transitory if 



H{x,y,t) 



H P (x,y), i<0 
H F (x,y), t>r 



If the transition time r goes to infinity, adiabatic theory [2^] may apply to such a 
system. To consider this limit, it is helpful to reformulate the equations by introducing 
a "slow time" variable 

A = et, 

where e = t" 1 is to be thought of as small. Instead of thinking of H as a function 
of time, through the transition function s, it is convenient to explicitly write it as a 
function of s, 

H(x,y,s(X)) = H(x,y,t), 

where s is a transition function with transition time 1. Now the past and future vector 
fields are generated by Hp = H(x,y,0) and Hp = H(x,y,l), respectively, and the 
vector field for (x, y, A) e M x 1 becomes 

d - 
x = —H(x,y,s(X)), 



V = - 
\ = e. 



d ~ 
—H{x,y,s(X)), 



(2.12) 
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The incompressible fluid flow (1.8) is a transitory Hamiltonian system in this sense 
with H(x,y,s) = -ip(x,y,t). 



For any finite e, the transition map (1.4 1 is now to be thought of as the time - 
map 

T(x,y) = <pi (x,y). 



The frozen time case corresponds to (]2. 12 1 with e — 0; alternatively it can be viewed as 



a family of autonomous systems on the phase space M with Hamiltonians H (x, y, s), 
aG [0,1]. 



Adiabatic theory shows that, in certain cases, orbits of (2.12) with e< 1 evolve 
so as to remain on orbits of the frozen time system with fixed loop action, where the 
loop action for a closed curve C is 



J(C) = d> ydx, (2.13) 

i.e. the area enclosed. Suppose there is a region R C M in which every orbit of 
the frozen system is periodic and that there is a smooth family 7,,. C R, s E [0, 1] of 
periodic orbits of the frozen time systems with fixed action 

J (7 S ) = J . 

Thus 70 is a periodic orbit of P, and 71 of F and J(jo) — J("fi)- According to 
adiabatic theory, if (x,y) € 70 and each of the 7 S has a bounded period, then 

d(T(x,y),-yi)^0 as e ->• 0, (2.14) 

where d(z, $7) = inf^ S Q \\ z ~ CI! is the standard distance from a point z to a set Q. 
The point is that T approximately maps periodic orbits of P to periodic orbits of F 
with the same action, T"(7o) ~ 71, when £«1. 

Adiabatic invariance breaks down if the frequencies of the periodic orbits in the 
family 7,, are not bounded away from zero. This occurs, for example, when 7^ ap- 
proaches or crosses a separatrix of the frozen time system for some s. The resulting 
jumps in the action due to separatrix crossing were computed by [51 145) . 

The flux in the adiabatic limit was studied by Kaper and Wiggins [5S] under 
the assumption that the frozen systems H (x, y, s) have a smooth, compact family of 
saddle equilibria, p(s), with homoclinic loops F p f s ) ^ W u (p(s)) D W s (p(s)). These 
authors state that the normally hyperbolic invariant manifold in the extended phase 
space, 

A = {(p( s (Xj),\):\eR}, 

continues to a nearby normally hyperbolic invariant manifold, A £ , when e is sufficiently 



small. If the Melnikov function for (2.12) has a nondegenerate zero, then the stable 



and unstable manifolds of A £ intersect transversely for small e defining a lobe R(s). 
Kaper and Wiggins show that the lobe area in the adiabatic limit becomes 



where 



lim Area(i?(e)) = J max - Jmin, ( 2 -15) 



J max — max J(T p i s )), J m in = mill J(T p ( s \) 

«e[o,i] «e[o,i] 
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are the maximum and minimum of the areas contained in the frozen time homoclinic 
loops. The implication is that the region that is "swept" by the separatrix is filled by 
the lobe, as was argued in [T5] . 



We will use ( 2.15 ) when possible in our examples below to discuss the limit r — > oo; 
however, the assumption that the frozen time systems have a normally hyperbolic 
manifold of equilibria can be easily violated for a transitory system. 

3. Examples. We will now consider several examples of transitory systems and 
use the results of 32] to quantify transport between coherent structures. For all exam- 



ples we use the cubic transition function s(t) — siit/r) from (A.l), scaled so that r 
is the transition time. Points along invariant manifolds are advected using a Runge- 
Kutta (5,4) Dormand-Prince pair. We represent the initial manifolds by a collection 
of equally-spaced points and add points using an adaptive interpolation method simi- 
lar to that of Hobson [55] when neighboring trajectories separate beyond a prescribed 
threshold. Similarly, trajectories are removed if the spacing decreases below a smaller 
threshold. We find the heteroclinic orbits from the advected manifolds using bisection 
to locate intersections of the stable and unstable invariant manifolds. 

3.1. Rotating Double Gyre. The motion of a passive scalar in a two-dimensional 
incompressible fluid with the oft-studied double-gyre configuration is depicted in the 



left pane of Fig. 3.1 This configuration has been observed in both geophysical flows 
[461 II 1 j and experimental investigations of laminar mixing in cavity flows |101 132] . 
Here we consider a transitory flow that corresponds to a rotation of the two gyres by 
| about (s, §)• It is defined by the stream function 

i>(x,y,t) = (l-s{t))i> P + s{t)Tp F , 

ipp(x,y) = 8in(2irx)sin(ny), (3-1) 

ipF(x,y) = sin(7rx) sin(27ry), 



and is Hamiltonian, with H = —ip and equations of motion (1.8). Such a rotating 



regime could arise in a geophysical setting if the prevailing direction of the jet separat- 



ing the two gyres changed over a finite time interval. In terms of a cavity flow, (3.1 1 
would model a regime in which a flow driven by upward movement of the left and 
right walls transitions smoothly to a flow driven by rightward movement of the top 



and bottom walls (cf. Fig. 3.1). Taking a slightly different perspective, the transition 
in (3.1) could also be effected by modulating the flow boundary over the transition 
interval. Such a changing boundary is used as a mixing mechanism in closed pipe 
flows, as investigated in |27|l3]. and is the driving force behind in-pipe mixing devices 
such as the Kenics® static mixer. 



The dynamics of (3.1 1 preserves the boundaries of the square M = [0, 1] x [0, 1] 
and each of the corners of M is an equilibrium of V. The corners (0,0) and (1, 1), 
are hyperbolic saddles for the full vector field V; however, though the corners (1,0) 
and (0, 1) are saddles for both P and F, they are not hyperbolic under the full field 
V. For example, the unstable manifold of (0, 1) is a subset of its stable manifold; the 
slices of these manifolds at t — are 

M/ U (0, 1) = {(x, 1) : < x < §} C 1F S (0, 1) = {(x, 1) : < x < 1}. 

Similarly the stable manifold of (1,0) is a subset of its unstable manifold; the slices 
at t = t are 

1^(1,0) = {(l,y) : < y < 1} C l^ r u (l,0) = {(l,y) : < y < 1}. 
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Fig. 3.1. Orbits of the stream functions ipp (left) and iftp (right) for \3. i[ ) . The unstable 
manifolds for the saddles of tpp are shown in red and the stable manifolds of the saddles of ipp are 
shown in blue. 



Thus, though (0, 1) and (1, 0) are both forward and backward hyperbolic, in the sense 
of Def. |1.2| neither is a hyperbolic orbit of V. 

The past vector field also has two saddle equilibria at po = (|,0) and pi = (i , 1), 
and the future vector field has saddles at /o = (0, |) and f\ = (1, ^). Under ( 3. lh the 
orbits of these points remain on the invariant boundaries of M and, as we shall see, 
these orbits play crucial roles in the delineation of Lagrangian coherent structures for 
this system and the quantification of the flux between them. 



The natural coherent structures for the past vector field of (3.1) are the left and 
right gyres separated by the unstable manifold 

W = {(i,y):0<y<l} (3.2) 

U. Moreover, for any t, the 
< x < 1}, and consequently 
7(pi,0) is a hyperbolic orbit of V. For the future vector field the top and bottom 
gyres are coherent structures with separatrix given by the stable manifold 



of pi, see Fig.[3Tj Note that for any t < 0, W t u (70i, 0)) 
stable manifold of this orbit is W^( r y(pi, 0)) = {(x, 1 



S={(xA) :0<x< 1} 



(3.3) 



of the future hyperbolic point f\. Note that W t s (7(/i, r)) = S for t > t, and 
W?( 1 (f 1 ,T)) = {(l,y):Q<y<l}. 

Transport between the two pairs of gyres is completely determined by the image 
of the past separatrix, T(U), or equivalently, the preimage of the future separatrix, 
T _1 (5). A movie showing the evolution of the manifolds starting with U and T~ 1 (S) 
at t = and ending with T(U) and S at t = r = 0.7 is at [LINK MOVIE "TDGAd- 
vection.mov" HERE]. However, to use the formulas of ^2l we only need to know the 
orbits heteroclinic from 7(pi,0) to 7(/i,t), and at t = r these consist of points on 
T(U) H S. At least one such heteroclinic orbit is guaranteed since dM is invariant 
under (3.1 1; that is, the segment T(U) still connects the top to the bottom and thus 
must cross S an odd number of times. When r = the transition map is the iden- 
tity and the only intersection is T(14) OS = hi = (|, 5) (see the top left pane of 
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Fig. 3.2. T(U) for transition times r = 0.0 (a), 0.4 (b), 0.5314 (c) and 0.8 (d) for the model 
\3. 1\ . Trajectories that begin in the left gyre are colored red, those that begin in the right gyre are 
colored blue, and the dividing curve between these is T(U). The dark blue region, labeled A r t, is 
the portion of the right gyre that ends in the top gyre at time r, and the dark red region, labeled 
An,, * s the portion of the left gyre that ends in the bottom gyre. The heteroclinic points are labeled 
{hi,h2,h:}} = T(U) n S. A movie showing the variation of T(U) as r varies from to 3.0 is at 
[LINK MOVIE "TDGManifolds.mov" HERE]. 



Fig. 3.2). This intersection persists as r increases, simply moving to the right along 



S and finally limiting on /j as r — > oo. For the cubic transition function, we observe 
that hi is the only intersection providing r < 0.5314, at which point a new pair of 
heteroclinic points, /12, ^3, are created in a saddle- node bifurcation, as shown in the 
bottom left pane of Fig. 3.2 The next heteroclinic bifurcation occurs at r sa 3.6908 
and the creation of heteroclinic orbits accelerates as r increases; indeed for r = 5.0 
there are 19 such heteroclinic orbits and the number appears to grow without bound 
as r — > 00. 



For the model (3.1) it is easy to obtain the manifold structure at t = from 



that at t = r, as shown in Fig. 3.2 because the system has a time-reversal symmetry 
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whenever the transition function obeys the relation 

s(r-t) = l-a(t). (3.4) 

Note that the cubic function we use and all of the polynomial transition functions 
given in App. \K\ have this property. In this case, the dynamics is reversed by the 
involution R : M x R ->• M x R defined by 

R(x,y,t) = (y,x,T-t). (3.5) 

It is easy to see that ip is invariant under this transformation and that the vector field 
is reversed by i?* , the push-forward ( |B.3[ ) of R: 

DRV(y,x,T-t) = -V(x,y,t). 

Consequently R inverts the transition map, T _1 = R o T o j?, and since R(S) = U, 

T- 1 {S)=R{T{U)). 

Consequently, the phase portraits of U and T~ 1 (S) at t = can be obtained from 
those in Fig. |3.2| at t — r by reflection about y = x and exchanging U and S. 

Denoting the left and right gyres of the past vector field by I and r and the top 
and bottom gyres of the future vector field by t and b respectively, there are four 
fluxes of interest, Aij, corresponding to the trajectories starting in gyre i £ {l,r} at 



time and ending in gyre j € {t, b} at time r (see Fig. 3.2). For example, A rt is the 
area of the region that is to the right of U for t < and above S for t > r. Thus, 
at t = t these regions are bounded by segments of T(U) and S. They can consist of 
multiple disjoint lobes when there are additional heteroclinic points, as in pane (d) of 
Fig. |3.2| in which there are two lobes for A rt . 



Since (3.1| is incompressible and dM is invariant, 

A n + A w = An + A rb = §, 

the total area of a single gyre in the past and future vector fields. Moreover, since the 
top and bottom gyres are filled completely by the images of the left and right gyres, 

Ait + A rt = An, + A r b = 2- 

Consequently, knowledge of one of these four areas uniquely determines the remaining 
three. 

To calculate A r t for a given transition time r we must integrate the form u> over 
the the dark blue region in Fig. |3.2| or equivalently, integrate the form —v over its 
boundary. We first consider the case of only one heteroclinic point h\ = T(U) n S. 
In this case, A r t is comprised of a single lobe and, given the points T(p\) = (x p , 1) 
and hi = (xh, \), the integrals along the top, right, and bottom edges of this lobe are 
trivial. Then, 

Kt = ^(l + Xh)-x p - / v, (3.6) 

where T(U^) is the oriented segment of T(U) from 7 r (pi,0) = T(p{) to 7 T (/ii,r 



h\. According to (2.4 1, the last term in (3.6) is simply the backward action difference 



AA T (T(pi),hi) between the orbits of h± and T(jp{). Thus, to evaluate (3.6) we 
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must compute T(jpx) and hi, and finally integrate the Lagrangian along the backward 
asymptotic orbits of hi and T(pi) from — oo to r to compute this action difference. 

Computation of T{jpi) is straightforward and is accomplished by numerically in- 
tegrating the vector field V over the transition interval [0, r] with initial condition 
pi. Computing hi is slightly more involved and requires a root-finding algorithm to 
determine the intersection of T(U) with S at t = t. We begin with points equally 
spaced along U at t — and numerically integrate V over [0, r] to compute the image 
under the transition map T of each point. As the manifold stretches during advection, 
we adaptively refine it to maintain the initial resolution. That is, at the first time 
step in which two neighboring trajectories diverge beyond a prescribed separation 
tolerance, we initialize a new trajectory at their midpoint and continue to track it 
over the remainder of the interval. We also remove points in much the same manner 
in areas where the manifold is contracting, helping to speed up computation. Upon 
obtaining T(U), we bracket the intersection with S and use a bisection routine to 
determine hi to the desired accuracy, initializing new trajectories where necessary. 
It should be noted that since unstable manifolds attract orbits in forward time, this 
adaptive refinement is a stable process and we incur minimal numerical error in the 
resulting manifolds T(U) [55]. In fact, results for T(U) obtained by refining the initial 
point spacing at t = and by adaptively refining during integration were virtually 
indistinguishable, with the adaptive computation being an order of magnitude faster 



in some cases. Finally, the Lagrangian (2.2 1, which for this system is 



L{x, y, t) = yx + i/j(x, y, t) = (1 - s{t))L P {x, y) + s(i)L F (x, y) 
Lp(x,y) = sin(27rx) [sin(7ry) — ny cos(7rj/)] , 
Lp(x,y) — sin(7ra) [ sin(27ry) — 27rycos(27ry)] , 



(3.7) 



is integrated along the computed orbits 7(pi,0) and "/(hi,r) using Simpson's rule 

When there are additional heteroclinic points (e.g., hi and ft. 3 in Fig. 3.2 
(d)), they can be computed in precisely the same way as hi, described above, 
the case of three heteroclinics the total flux is A rt = A\ t + A^, where the two terms 
on the right-hand side represent the areas of the two disjoint lobes. The area of the 



pane 
For 



larger lobe A\ t is calculated using (3.6) while the smaller lobe is similar to that shown 
in Fig. 



2.1 and hence we calculate its area according to (2.7) 



/oo 
[L{h 2 (t),t)-L(h 3 {t),t)]dt, 
-oo 



which is positive by the counterclockwise orientation of the segments U h l and S h 



Computation of the last term in (3.6) is grea tly simplified using (2.9), since the 
system is transitory, and by noting that for (3.7) Lp{h,y) = for any y. Similarly, 
(2.11 ) can be used to simplify the computation of A^f. 

Results for the computation of A rt for transition times ranging from to 3.69 are 
summarized in Fig. |3.3| Note the increase in the rate of change of flux at r w 0.531 



corresponding to the emergence of the second lobe of area A^f. At r « 3.69, a new 
pair of heteroclinic points h^ and h$ is created and their corresponding manifolds 
delineate a new lobe of trajectories initially to the left of U for t < 0. Indeed, each 
new heteroclinic bifurcation as r increases creates a new lobe that alternately adds 
area to A rt or to A\ t . 
As r 



oo some of the orbits of (3.1) can be described by the adiabatic theory 



outlined in £2.2 For example, each periodic orbit in the neighborhood of the elliptic 
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Fig. 3.3. Plot of the right-to-top flux, A r t, as a function of transition time. Note that a second 
lobe of area A^ emerges at r cr i t ~ 0.531. Its effect is manifested by the departure of the solid line 
from the dashed one for A^, t . 



equilibrium ol the left gyre of P continues to a periodic orbit of the frozen-time 
system — setting ip(x,y, s(t)) = tjj(x,y,t) — with fixed loop action. As s grows from 



to 1 in (3.1 ) the left gyre rotates by § , so these orbits evolve continuously to periodic 
orbits of F enclosing the elliptic equilibrium of the bottom gyre. If a family of periodic 
orbits of the frozen time system with fixed loop action remains a bounded distance 
away from the family of separatrices of ?p(x,y,s), then the period of each orbit in 
this family is bounded, and as r — > oo the adiabatic theory implies that the actual 
evolution will follow that of the frozen system. The implication is that when r ^> 1, 
T will approximately map periodic orbits of P in the left gyre to periodic orbits of F 
in the bottom gyre with the same action (and similarly for the right and top gyres). 
An indication of the approach to adiabaticity is displayed in Fig. |3.4| which shows 
four elliptic orbits 7q, i = 1, ..., 4 of P (left pane) and their images under the transition 
map T for two values of r (middle and right panes). The dashed curves represent 
orbits j\ of F having the same loop actions as the initial orbits. When r is small, 
as in the middle pane, each of the images T(jq) differs visibly from r )\. Conversely, 
when r is moderately large, as in the rightmost pane, T maps the innermost three 
7g virtually on top of the corresponding r y\. The outermost orbit, 7q, (red curve) 
maps to a loop with tendrils far from jf as its period is not sufficiently small for 
adiabaticity to pertain at this value of r. Violation of adiabaticity could also be 
seen in its most extravagant form by T(U) itself, which, as noted above, necessarily 
crosses the square from top to bottom and intersects the line S infinitely many times 
as r — > oo. Nevertheless, since increasingly many orbits of the left gyre map to their 
counterparts in the lower gyre as r increases, adiabatic theory implies that 

An, = A rt — > 0.5 as r — > oo, 

as Fig. |3.3| seems to suggest. 

While this result is clearly demonstrated by the numerical evidence shown in 
Fig. |3.4| and the corresponding movie linked in its caption, the theory of Kaper and 
Wiggins jUj for the flux in the adiabatic limit, outlined in { 2.2 above, does not apply 



directly to the double-gyre model (3.1). There is no family of homoclinic loops for 
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i/)(x, y, s(A)), though, as Kaper and Wiggins themselves point out, the theory could be 
straightforwardly extended for a family of heteroclinic cycles, as would be appropriate 
for the double-gyre. When s < |, the "left" gyre of ip(x 7 y,s(\)) is bounded by 
separatrices connecting four hyperbolic equilibria: (0,0), (0, 1), and a saddle on each 
of the upper and lower boundaries. This family of separatrices looses hyperbolicity 
at s = 27 when the point (0,1) is no longer a hyperbolic equilibrium of the frozen 
system and the separatrix becomes a triangle. Though the p oint (0 , 1) again becomes 
hyperbolic when s > 
applied. 



it does not appear that the result (2.15) can be rigorously 




Fig. 3.4. Illustration of the approach to adiabatic invariance of periodic orbits for (3. 7| ). The 
solid curves depict the images of four orbits of P under the transition map T for r = (left pane), 
t = 0.3 (middle pane), and r = 2.5 (right pane). The dashed curves represent orbits of the future 
vector field with the same loop actions as the initial orbits. The solid black curve depicts the orbit 
of the left elliptic equilibrium of P over < t < r. A movie showing the emergence of adiabatic 
behavior as r increases is [LINK MOVIE "TDGAdiabatic.mov" HERE] 



Finally, we compare our mode of analysis with a technique commonly used for 
analyzing time-dependent flows: the finite-time Lyapunov exponent (FTLE). We will 
comment only briefly here on the similarities and differences between these two meth- 
ods in the context of identifying heteroclinic trajectories and computing lobe areas for 
(3.1). A more thorough explanation of the use of FTLE for approximating invariant 
manifolds is given in [52] [18] . 

The backward-time FTLE field for (3.1) with transition time r = 0.8, is shown 



in the left pane of Fig. |3.5| The "ridges" of the FTLE field approximate the unstable 



boundaries of LCS; the red regions in the figure. Comparing Fig. 3.5 with pane (d) of 
Fig. |3.2[ one can see that the most prominent ridge of the FTLE field corresponds to 
the curve T(U). However, the global picture given by the FTLE field is complicated 
by secondary ridges near the main ridge. These secondary ridges are common (see for 
example [5211181 17]) and, while it is unclear whether they are numerical anomalies or 
they offer physical insight into the stretching of nearby trajectories over the time-scale 
used for computation, they complicate the numerical extraction of the most prominent 
ridge and the identification of any heteroclinic orbits. 

As noted in [52] , the ridges of the FTLE field become "more Lagrangian" as the 
integration time grows. Since the transitory system (3.1) has a trajectory ipt,o(Pi) 
that is truly backward hyperbolic (i.e. it does not lose its backward hyperbolicity 
for any time t € K), as the integration time increases we should expect the most 
prominent ridge of the backward-time FTLE to become increasingly aligned with the 
unstable manifold T{U) of T(p\) at time r. We do indeed observe this; however, 
the secondary ridges also become more pronounced, making extraction of the main 
ridge increasingly difficult. This places a practical upper limit on the length of the 
approximate invariant manifold that can be computed using FTLE calculations, as 
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Fig. 3.5. Backward time FTLE field for \3.1\ at the transition time r = 0.8 using a backward 
integration time of 1.2 and a 1500 X 1500 grid. The left pane shows contours of the FTLE value 
from blue (smallest) to red (largest), and the right pane shows the ridge extracted by keeping only 
those values within 30% of the maximum. This ridge gives an approximation to T(U) . 



the numerical extraction of the "main" ridge becomes infeasible for large integration 
times. 

Several numerical methods for efficiently extracting the appropriate ridges from 
the FTLE field have been proposed J49J, HTJ, [33] ; however, in practice, the most promi- 
nent ridges are typically extracted by simply filtering out all values below a prescribed 
threshold. Such a filter with the threshold set at 70% of the maximum FTLE value 
is shown in the right pane of Fig. 3.5 For the chosen integration time, the secondary 



ridges are so close in height to the main ridge that they can not be removed by this 
simple height filter. That is, as the filtering threshold is increased, gaps appear in 
the main ridge before all the secondary ridges have disappeared. Thus, while the 
main FTLE ridge qualitatively agrees with the true unstable manifold T(jU), an ac- 
curate computation of flux between coherent structures is difficult to obtain with this 
method. 

3.2. Resonant Accelerator. As a second example, we consider a system that 
serves as a highly simplified model of a particle accelerator [T3]. Here the coherent 
structures are "resonances" that result in the trapping of particles in an accelerating 
potential well, and the goal is to determine the phase space region that represents 
stable acceleration. In our model, this corresponds to orbits that begin within a 
stationary, past resonance and ultimately end in a moving, future resonance at t = r. 

The basic model is given by a Hamiltonian of the form 



H(q,p,t) = ^p 2 



V(q-9(t)). 



We assume that the potential well is initially stationary, then accelerates, and even- 
tually reaches a constant velocity so that the phase 9(t) obeys 



6{t) = 





LUt + (j) 



t<0 

t> T 



(3.8) 



While this system is not transitory in the sense of Def. |1.1| (note the time dependence 
of the potential function V), we can convert it to one that is with the canonical 
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transformation 

(q,p,H) k> (Q,V,H) = (q-e(t),p,H-p6(t)). 



Note from_J3.8) that the time derivative of the phase is proportional to a transition 



function (1.6 1, namely 9{t) — u>s(t). Reverting to the original variable names gives 



the new Hamiltonian 

H^p^^^-us^p + Viq), (3.9) 

and so the past and future systems are autonomous: 

H( q ,p,t) = ( ^ P \= |f + V #^ 12 I <0 • (3.10) 

w.i'. i \H F (q,p)= \{p-uy + V{q)-\u 2 , t>r 

Specifically, for our model we take 

V(q) = — fccos(2-7rg), 

so both autonomous limits are equivalent to the pendulum, with the resonance cen- 
tered around p = for t < and p = ui for t > r. Without loss of generality we can 
scale variables to set u = 1, leaving two parameters: the transition time r and the 
potential energy amplitude k. All examples and parameter values below correspond 
to this scaled system. 

Contours of Hp and Hp are shown in Fig. |3.6[ and since oj = 1, the transition 
Up — > ifp corresponds to a unit vertical translation of the past vector field. For 
each s, the frozen-time Hamiltonian has saddles at (q,p) = (±i,s) with stable and 
unstable manifolds 

p±(q, s) = ±y/2k(l + cos(2nq)) + s, (3.11) 

that define separatrices bounding a resonance. Of course, by periodicity the two 
saddles can be identified, so that the manifolds actually correspond to homoclinic 
loops on the cylinder M = S X K. The past resonance corresponds to s — and its 
separatrices are denoted U± as these are slices of the unstable manifolds of the saddle 
for the full vector field when t < 0. The future resonance corresponds to s = 1 and 
its separatrices are denoted S± , as these are slices of stable manifolds of the saddle 
for the full vector field when t > r. The width of each resonance of the frozen time 
system is 

w=p + (0,s)-p-(0,s) =4\/fc, (3.12) 

which for this simple model is independent of s. Note that the past and future 
resonance zones "overlap" when w > 1, implying that k > A. 

Let Ai be the area of the region that begins inside the past resonance at time 
t = and ends outside the future resonance at time t = r, with corresponding 
notations A^ , A oi and ^4 00 for the other beginning and ending configurations. We are 
principally interested in calculating the fraction of accelerated phase space area, 

A.. 
Rac 



The images of the manifolds IA± under the transition map for two values of r are shown 



in Fig. 3.7 Here Ai is the area of the region that is inside T{U±) and outside S±. 
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Fig. 3.6. Contours of the past (left) and future (right) Hamiltonians f 3. 1 0[ ) for k = 0.4. The 
unstable manifolds for the saddles of the past vector field are shown in red and the stable manifolds 
for the saddles of the future vector field are shown in blue. 
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Fig. 3.7. Image of the unstable manifolds of the past resonance for \3.Sty at time r for k = 0.4 
with transition times r = 1.0 (left) and r = 3.0 (right). The blue regions correspond to trajectories 
that begin inside the past resonance and the red regions to those that begin outside. The dark blue 
region, labeled An corresponds to those particles that remain trapped in the accelerated potential well 
and the light blue region, Ai to those that are left behind. 



This region, light blue in the figure, appears disconnected; however, on the cylinder 
it is formed from a single connected set. The region An is dark blue in the figure and 
corresponds to particles that remain trapped within the resonance for t > r. 

In Fig. 3.7 there are two heteroclinic points, {hx,fi2} € T(U- 
such heteroclinic orbits do not always exist. In particular, for r ; 



n S- ; however, 
the transition 



map is the identity and so there are heteroclinic points only when the past and future 
resonances overlap, that is, when k > k cr u(0) = A. More generally, the critical k for 
a heteroclinic bifurcation can be computed as for the double gyre model (3.1). The 
resulting curve of bifurcations, k cr it(j)i is shown in Fig. 3.8 and a phase portrait at 
the bifurcation point for t = r = 1 is shown in the right pane of this figure. When 
k > kcritij), there is a pair of heteroclinic orbits, and — unlike the double gyre — there 
appear to be no additional heteroclinic bifurcations as r grows. Given the orbits of h\ 
and /i 2 , the area A^ can be computed according to the simplified lobe area formula 
(2.11 ), since the system is indeed transitory. The resulting ratio R aC c as a function of 



both k and r is shown in Fig. 3.9 
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It is interesting that we never observe intersections between T(U±) and <S+. In- 
deed, it is easy to see that there is no flux out the top of the instantaneous separatrix. 
To see this, let 

F S ep{ s ) = n ~ 2 S ' 

be the energy of the frozen separatrix and define 

F{q,p, s) = H(q,p, s) - E sep (s). (3.13) 

Note that F < inside the separatrix and F > outside of it. Differentiation along 
the Hamiltonian vector field gives 

F = —s(p — s). 

Since p > s for any point on the separatrix p + (q,s) and since we assume s(t) > 0, 
then F < on the upper separatrix. Thus the vector field never permits a trajectory 
to cross this separatrix from below, and a transverse intersection of T{U±) with <S + 
is forbidden. 




Tradition Time (t) 




Fig. 3.8. Curve of heteroclinic bifurcations for \3. 9[ ) (left pane) and the bifurcation point for 
t = t = 1.0 with k = k cr it ~ 0.057 (right pane). 



The unstable manifold in the right pane of Fig. 3.7 can be compared with the 
ridges of the backward-time FTLE field in Fig. |3.10[ As in the double-gyre example, 
the qualitative agreement is quite good, even though secondary ridges do exist. A 
simple threshold filtering of the FTLE field, shown in the right pane of the figure, 
extracts the most prominent ridge; however, this set is not a curve on the scale of 
the computational grid, especially in the interior and near the boundary of the future 
resonance. Once again, the secondary ridges are so similar in height to the main 
ridge that they can not be removed by a height filter. This ambiguity as to the true 
location of the unstable manifold makes it hard to identify the heteroclinic points 
and accurately compute the desired flux from the manifolds obtained from the FTLE 
field. 

Adiabatic theory applies to the resonant accelerator system once it is written in 



the form (2.12). There are two topologically distinct loops in the frozen-time phase 



space corresponding to trapped (oscillatory) and to untrapped (rotational) trajecto- 
ries. For the trapped case the loop action is the area enclosed, but for the untrapped 
orbits the loop action is the area contained between the graph of the trajectory and 
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Fig. 3.9. Fraction of the past resonance that is trapped in the future resonance for \3.9\ as a 
function of parameters k and r. 




-0.6 -0.4 -0.3 -11.2 -0.1 » 0.1 11.2 11.3 11.4 0.8 



-0.5 -0.4 -0.3 -11.2 -0.1 II D.l 0.2 0.3 11.4 0.5 



Fig. 3.10. (left) Backward-time FTLE field for 1(3.9]) with k = 0.4, r = 3.0 and integration 
time t = —6. The black line denotes the resonance of the future vector field F. (right) The main 
ridge extracted from FTLE field by keeping only those values within 45% of the maximum. 



the g-axis. Since k is constant in our model, the area of the trapped region of the 
frozen system is independent of s, and thus for each trapped orbit of P, there is a 
family, 7 S , of periodic orbits of the frozen systems with the same action. Each trapped 



orbit that is a bounded distance inside the separatrices (3.11) has a period that is 
bounded; consequently, these orbits will be adiabatic in the limit r — > oo. For the 
untrapped orbits, this is no longer always true. As s varies, the upper separatrix of 
the frozen system "sweeps through" the region bounded by the separatrix IA + , with 
action J (14+) = -y/2k, and the separatrix <S + , with action J{S+) = J(U+) + 1. Every 
rotational family of orbits of the frozen system within this separatrix- swept region, 
namely those with actions 

Je[J(U+),J(S + )}, 

necessarily crosses the separatrix p+(q, s) for some s £ [0, 1]. Since this implies the 
period is unbounded, adiabatic theory does not apply to these orbits. 
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By contrast the rotational orbits with J < J{U-) or J > J{S+) remain bounded 
away from the frozen separatrices for all s, and thus are adiabatic in the limit r — > oo. 
This is supported by the computations in Fig. |3.1l] and the corresponding movie linked 
in its caption. The left pane of the figure shows eight orbits of P, and the middle 
and right panes show the images of these under T for two values of r. When t = 10, 
each of the orbits — except for the green orbit in the separatrix swept region — has an 
image under T that is very close to an orbit of F with the same loop action. 




0.2 -0.1 0. 



Fig. 3.11. Illustration of adiabatic invariance of the loop actions of periodic orbits for 13. 9[ ) . 
The left pane shows eight orbits 7q of the past vector field, and the middle and right panes depict 
T(7q) for t = 1.5 and 10. The blue curve is the boundary S± of the future resonance and the dotted 
curves indicate orbits of F with the same actions as the orbits ■Jq . The black curve shows the orbit 
of the elliptic equilibrium of P. A movie showing the variation of these images with r is at [LINK 
MOVIE "AccAdiabatic.mov" HERE] 



The frozen-time accelerator model has a family of hyperbolic saddles, and each 
saddle has a pair of homoclinic loops. Thus, the theory of [25] described in £2.2 
applies. Since J max — J m in for this model, the area of the lobe, Ai , limits to zero, 
and so this theory predicts that 

lim R ar r = 1. 



That is to say, the region of area An limits to the future resonance, as Fig. 3.7 and 
the accompanying movie seem to suggest. In terms of the physical model this implies 
that, provided the acceleration is "slow enough," almost all particles beginning in the 
resonance at t = will be stably accelerated. 

4. Conclusions. While techniques involving finite-time Lyapunov exponents 
and distinguished hyperbolic trajectories have recently been developed for the iden- 
tification and extraction of coherent structures in time-dependent systems, they have 
been used only selectively to give quantitative descriptions of the finite-time flux be- 
tween such structures. One reason for this is that the "ridges" of the FTLE field 
that represent approximate invariant manifolds are often difficult to extract, making 
precise measurements of flux challenging. 

Here we have considered a special class of two-dimensional nonautonomous sys- 
tems that exhibit time-dependent behavior only on a compact interval, and have 
extensively used the concepts of backward and forward hyperbolicity for these tran- 
sitory systems. The special structure of these systems, leads to a simple a method 
for the numerical computation of flux between Lagrangian coherent structures in the 
Hamiltonian case. Our method relies primarily on knowledge of heteroclinic orbits 
and their associated invariant manifolds that bound lobes within the extended phase 
space. Thus, our computations of flux require very little Lagrangian information rel- 
ative to computations involving FTLE or distinguished hyperbolic trajectories. In 
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particular, our adaptive computation of T(U) allowed for an order of magnitude re- 
duction in the number of particle advections required for a computation of the FTLE 
field at similar resolution. 

An important extension to the theory presented here arises in light of recent 
advances in the theory of finite-time manifolds [20l EU [12 [54] . These studies have 
shown that, given a system whose behavior is unknown outside an interval X = [£_ , t+] , 
the manifold structure of a "sufficiently slowly" moving orbit at some time t £ int(I) 
is "unique" up to an exponentially small correction term, provided t is "sufficiently 



far" from the endpoints of I. Since the formulas (2.4) and (2.6| depend only on 
heteroclinic orbits lying on the boundary of a lobe, they could be used directly to 
provide exponentially accurate approximations of lobe areas in such systems. 

It is not obvious if our technique can be applied more generally to nonautonomous 
systems or to systems defined by a discrete set of data; however, there are several 
logical extensions that we plan to address in future work. The first is the case of tran- 
sitory, symplectic maps, where the action formulas that we have developed should also 
apply. Similarly, since action formulas for flux have been recently developed in the n- 
dimensional volume preserving case 34J, transitory volume-preserving systems could 
also be treated. Finally, it is reasonable that if the time dependence is "episodic" in 
nature, each transition could at least approximately be treated by the same methods 
that we have used. 



Appendices 

Appendix A. Transition functions. 

We can obtain a C k transition function that is polynomial ont £ [0, 1], by requir- 
ing s(0) = 0, s(l) = 1 and ZPs(O) = D j s(l) = for j = 1 . . . k. On the interval [0, 1] 
these functions are 

s (t)=t 

Sl (t) =£ 2 (3-2f) 

s 2 (t) = i 3 (10 - 15t + 6i 2 ) (A.l) 

s 3 (t) = i 4 (35 - 84t + 70i 2 - 20i 3 ) 

Si (t) = i 5 (126 - 420t + 540t 2 - 315i 3 + 70t 4 ). 

Figure PTT] shows several of these functions for various values of k. It is not hard to 
see that in general these polynomials are given by 



r (») /'.», 



" (t) -(f5#X' (1 - 8)d " 

for t £ [0, 1], which is monotone. 

Appendix B. Forms and Lie Derivatives. Here we set out our notation, 
which follows [T]. We denote the set of k- forms on a manifold M by A fc (M), and 
the set of vector fields by V(M). If a £ A k (M) and Vi, V 2 , ■ ■ . V k £ V(M), then the 
pullback, /*, of a form a by a diffeomorphism / is defined by 

(f* a ) x (V 1 ,V 2 ,...,V k ) = a nx) (Df(x)V 1 (x),...,Df(x)V k (x)). (B.l) 

The pullback can be applied to a vector field V as well: 

(f*V)(x) = (Df(x))-W(f(x)). (B.2) 
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Fig. A.l. The transition functions sj,(t) for odd k up to 9. 



(B.3) 



The push-forward operator is defined as 

/* = (r 1 )*. 

The interior product of a with V is defined as the (k — l)-form 

i v a = a{V,-,. ..,■). (B.4) 

Suppose that ip t ,t is the (C 1 ) flow of a vector field V(x,t), so that ip to j (x) — x, and 
iiVt,t {x) = V(ip t ,t 
operator defined by 

d 



^ip ti t (x) — V(ipt,t (x),t). Then the Lie derivative with respect to V is the linear 



C v - = 



ds 



v s , t - 



where • is any tensor. In particular for a vector field X, 

£ V X = [V, X] = (V ■ V)X - (X ■ V) V, 



(B.5) 



(B.6) 



where [ , ] is the Lie bracket. The Lie derivative acting on differential forms obeys 
Cartan's homotopy formula 



£yQ = %v{da) + d(lyOt). 
Note that C behaves "naturally" with respect to the pullback: 

f*C v a = Cf*vf*a. 



(B.7) 



(B. 
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